/// Inequlities Solver


typedef enum {TREM = 3, TGEQ=2, TGE=1, TEQ=0, TINCONS=4} Constraint_Type;
#define DEBUG_FM_SOLVER

template <class Matrix_ptr, class Vector_ptr, class Types_ptr>
class FM_Solver {
    Matrix_ptr constraints_;	/// Pointer to the matrix of inequlities
    Vector_ptr constants_;	/// Pointer to the vector of constatnt
    Types_ptr  ctypes_;		/// Pointer to the vector of types
public:
    typedef std::set<int> Vars;
    FM_Solver();
    FM_Solver(Matrix_ptr m, Vector_ptr c, Types_ptr t);
    ProjectVar(int var_id);
    ProjectVars(Vars & var_ids);
private:
    /// Build projective matrix
    BuildProjective(size_t column);
    /// Clean up routines
    
    /// choose best to eliminate
    choice_inq(const Vars &vars)
    
};



Matrix_ptr FM_Solver::BuildProjective(size_t column) {

#ifdef DEBUG_FM_SOLVER
    cout<<"\nProblem::BuildProjective: column = "<<column;
    if (constraints_) cout<<"\nconstraints_ is valid"; else cout<<"\nconstraints_ is invalid";
#endif
    assert(column<constraints_->size2());
    Coefs varCoef;
    size_t size1 = constraints_->size1();
    
    for(size_t s1=0;s1<size1; s1++)
    	varCoef.push_back((*constraints_)(s1,column));
    
    PInts ref_pos; // index of positive members of varCoef
    PInts ref_neg; // .. negative ..
    PInts ref_nul; // .. 0 .. or other than GEQ, GE

    // buld indexes
    for (size_t i = 0; i<size1;i++) {
        if (varCoef[i]>0) {
            ref_pos.push_back(i);
        } else if (varCoef[i]<0) {
            ref_neg.push_back(i);
        } else {
            ref_nul.push_back(i);
        }
    }
    //if only negative of only positive exists we
    //just delete column in matrix A
    // Also We delete this colums anyway
    // constraints_->erase_line2(column);    /// smx1
    constraints_->reset_s2(column);     /// smx2 does not have erase_line2 method
    //if we have only one type of inequality
    //we mark it as REM in ctypes_
    // letter it will be deleted by remove_eliminated()
#ifdef DEBUG_FM_SOLVER
    debug("\npositives:", ref_pos);
    debug("\nnegativess:", ref_neg);
    debug("\nnuls:", ref_nul);
#endif
    if (ref_pos.empty() || ref_neg.empty()){
        std::for_each(
            ref_pos.begin(),
            ref_pos.end(),
            boost::bind(&Problem::reset_line,this,_1)
         );	
         std::for_each(
            ref_neg.begin(),
            ref_neg.end(),
            boost::bind(&Problem::reset_line,this,_1)
         );
        return SparseMatrix_aptr(); // empty
    }
    // reseting
    /// TODO: profile,and if too many ref_nuls
    /// move those ref_nul to other system
    /// and merge it to the main system after projection
    size_t new_size = ref_pos.size()*ref_neg.size()+ref_nul.size();
    debugf("\nprojective matrix(", new_size, ", ", size1, ") created.");
    //cout<<"\rmatrix("<< new_size<<", "<<size1<<") created."<<flush;
        
    if (new_size>=INT_PROBLEM_SIZE){
    	cerr<<"Maximum allowed internal problem size: "<<INT_PROBLEM_SIZE;
    	cerr<<"\nThe problem too complicated ("<<new_size<<", "<<size1<<")\nAborted."<<endl<<flush;
    	exit(1);
    }
    /**
    For sake of perfomance guessing of new number of
    non empty elements for the projective matrix is useful to
    implement.
   */
    SparseMatrix_aptr projective(new SparseMatrix(new_size, size1));
    ConstraintsType_aptr newctypes_(new ConstraintsType(new_size));	
    //buildin C for positive and negative
    // assuming that ref_pos and ref_neg at least has size 1
    int CurRow=0;
    for(size_t pp=0; pp < ref_pos.size(); pp++ ) {
        for(size_t pn = 0; pn < ref_neg.size(); pn++) {
            (*projective)(CurRow,ref_pos[pp]) = -varCoef[ref_neg[pn]];
            (*projective)(CurRow,ref_neg[pn]) = varCoef[ref_pos[pp]];
            // note. there are no equalities here
            if ((*ctypes_)[ref_pos[pp]] == (*ctypes_)[ref_neg[pn]]){
            	(*newctypes_)[CurRow] = (*ctypes_)[ref_pos[pp]];
            } else {
            	if ((*ctypes_)[ref_pos[pp]].type == TGE)
            		(*newctypes_)[CurRow] = (*ctypes_)[ref_pos[pp]];  // > * >= substituted by >
            	else {
#ifndef NDEBUG	
            	     if ((*ctypes_)[ref_pos[pp]].type != TGEQ){
            	     	cout<<"\nExpecting TREQ but got:"
            	     	    <<((*ctypes_)[ref_pos[pp]].type)<<endl<<flush;
            	     }
#endif
            	     assert((*ctypes_)[ref_pos[pp]].type == TGEQ);

            	    (*newctypes_)[CurRow] = (*ctypes_)[ref_neg[pn]];
            	}
            }
            CurRow++;
        }
    }
    // bilding C for 0 by seting it to 1
    // assuming size of ref_nul could be 0
    
    for(size_t p0 = 0; p0 < ref_nul.size(); p0++) {
        (*projective)(CurRow, ref_nul[p0]) = 1;
        (*newctypes_)[CurRow] = (*ctypes_)[ref_nul[p0]];
        CurRow++;
    }
    ctypes_ = newctypes_;
#ifdef DEBUG_FM_SOLVER    
    cout<<"\nbuilded c\n"<<*projective<<flush;
#endif
    return projective;
};
 
 void FM_Solver::EliminateVariable(int var_id) {
    SparseMatrix_aptr  projective = BuildProjective(var_id);
    if (projective) {
#ifdef DEBUG_FM_SOLVER
        cout<<"\ninq  "<<*constraints_;
        cout<<"\nconst  "<<*constants_;
        cout<<"\nproj "<<*projective;
    	cout<<"\nEliminateVariable before mult s11 "<<constraints_->size1()<<
    	" s21 "<<constraints_->size2()<<" s12 "<<projective->size1()<<
    	" s22 "<<projective->size2()<<endl<<flush;
#endif
        (*constraints_) = (*projective) * (*constraints_);
        projective->prod(*constants_);  /// constraints_ = projective * constraints_ 
        debugf("\nproblem size(",constraints_->size1(), ", ", constraints_->size2(), ")");
        assert(constraints_->size1() == constants_->size());
	
    } else {
        //alredy eliminated
    }
    
    return;
};

int Problem::choice_inq(const Variables &vars) const{
    // For each varible of type Wildcard_Var
    // calculate n_pos * n_neg
    // if the calulated value is less then current one
    // save the index and the value

#ifdef DEBUG_PROBLEM
    cout<<"\nchoice_inq: vars.size: "<<vars.size();
#endif
    int index=-1;
    if (!EnableVaribleChoise)
    	return (vars.size()? 0: index);
    
    int value = std::numeric_limits<int>::max();

    for(Variables::const_iterator it= vars.begin(); it!=vars.end(); ++it){
#ifdef DEBUG_PROBLEM
    	cout<<"\nchoice_inq: variable: "<<(*it)->name(); 
    	cout<<"\nchoice_inq: var Type: "<<(*it)->GetType();
    	cout<<"\nchoice_inq: var   id: "<<(*it)->id();  
#endif
        //if((*it)->GetType() != Wildcard_Var) continue;
        assert((*it)->GetType() == Wildcard_Var);
    	int c_val = var_difficulty((*it)->id());
#ifdef DEBUG_PROBLEM
    	cout<<"\nchoice_inq: Difficulty "<<c_val;
#endif
    	if (c_val < value){
    	    index = it - vars.begin();
    	    value = c_val;
    	    if (value ==0 ) // cannot be simpler so exiting
    	    	break;  
    	}
    }
#ifdef DEBUG_PROBLEM
    cout<<"\nchoice_inq: out: index "<<index;
    cout<<"\nchoice_inq: out: coml  "<<value;
#endif
    return index;
}
